function CeShen()
%给定参数
Nt = 80;
p1 = [100, 100, 10, 10, 10, 10, 100, 100];
p2 = [10, 10, 100, 100, 50, 50, 50, 50];
p3 = [50, 50, 50, 50, 100, 100, 10, 10];
h1 = 1000;
h2 = [200, 5000, 200, 5000, 200, 5000, 200, 5000];
u2 = p2 ./ p1; u3 = p3 ./ p2;
v2 = h2 ./ h1;
i = sqrt(-1);
%进行计算
for ix = 1:8
    for  iy = 1: Nt
        B(ix,iy) = 10 ^ ((iy+10) / 10) / h1;
        C(ix,iy) = ( coth( 2*pi*sqrt(2) * exp( -i*pi/4) / B(ix,iy) + ...
            acoth( sqrt(u2(ix)) * coth( 2*pi*sqrt(2) * exp( -i*pi/4) / ...
            B(ix,iy) * v2(ix) / sqrt(u2(ix)) + acoth( sqrt( u3(ix)/u2(ix)))))))^2;
        D(ix,iy) = ( tanh( 2*pi*sqrt(2) * exp( -i*pi/4) / B(ix,iy) + ...
            atanh( sqrt(u2(ix)) * tanh( 2*pi*sqrt(2) * exp( -i*pi/4) / ...
            B(ix,iy) * v2(ix) / sqrt(u2(ix)) + atanh( sqrt( u3(ix)/u2(ix)))))))^2;
    end
    Re1(ix,:) = real(C(ix,:));          %实部
    Im1(ix,:) = imag(D(ix,:));          %虚部
    A(ix,:) = abs(C(ix,:));
end
%进行画图
figure(1)
loglog(B(1,:),A(1,:))
title("H型三层大地电磁测深响应曲线(中间薄)")
xlabel("lanbta1/h1"); ylabel("pw/p1")

figure(2)
loglog(B(2,:),A(2,:))
title("H型三层大地电磁测深响应曲线(中间厚)")
xlabel("lanbta1/h1"); ylabel("pw/p1")

figure(3)
loglog(B(3,:),A(3,:))
title("K型三层大地电磁测深响应曲线(中间薄)")
xlabel("lanbta1/h1"); ylabel("pw/p1")

figure(4)
loglog(B(4,:),A(4,:))
title("K型三层大地电磁测深响应曲线(中间厚)")
xlabel("lanbta1/h1"); ylabel("pw/p1")

figure(5)
loglog(B(5,:),A(5,:))
title("A型三层大地电磁测深响应曲线(中间薄)")
xlabel("lanbta1/h1"); ylabel("pw/p1")

figure(6)
loglog(B(6,:),Re1(6,:))
title("A型三层大地电磁测深响应曲线(中间厚)")
xlabel("lanbta1/h1"); ylabel("pw/p1")

figure(7)
loglog(B(7,:),Re1(7,:))
title("Q型三层大地电磁测深响应曲线(中间薄)")
xlabel("lanbta1/h1"); ylabel("pw/p1")

figure(8)
loglog(B(8,:),Re1(8,:))
title("Q型三层大地电磁测深响应曲线(中间厚)")
xlabel("lanbta1/h1"); ylabel("pw/p1")

end
